/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003 RiskMap srl
 Copyright (C) 2015 Peter Caspers

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include "integrals.hpp"
#include "utilities.hpp"
#include <ql/math/integrals/exponentialintegrals.hpp>
#include <ql/math/integrals/filonintegral.hpp>
#include <ql/math/integrals/segmentintegral.hpp>
#include <ql/math/integrals/simpsonintegral.hpp>
#include <ql/math/integrals/trapezoidintegral.hpp>
#include <ql/math/integrals/kronrodintegral.hpp>
#include <ql/math/integrals/gausslobattointegral.hpp>
#include <ql/math/integrals/discreteintegrals.hpp>
#include <ql/math/integrals/tanhsinhintegral.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <ql/math/interpolations/bilinearinterpolation.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/termstructures/volatility/abcd.hpp>
#include <ql/math/integrals/twodimensionalintegral.hpp>
#include <ql/experimental/math/piecewisefunction.hpp>
#include <ql/experimental/math/piecewiseintegral.hpp>

using namespace QuantLib;
using boost::unit_test_framework::test_suite;

namespace integrals_test {

    Real tolerance = 1.0e-6;

    template <class T>
    void testSingle(const T& I, const std::string& tag,
                    const ext::function<Real (Real)>& f,
                    Real xMin, Real xMax, Real expected) {
        Real calculated = I(f,xMin,xMax);
        if (std::fabs(calculated-expected) > integrals_test::tolerance) {
            BOOST_FAIL(std::setprecision(10)
                       << "integrating " << tag
                       << "    calculated: " << calculated
                       << "    expected:   " << expected);
        }
    }

    template <class T>
    void testSeveral(const T& I) {
        testSingle(I, "f(x) = 0", [](Real x) -> Real { return 0.0; }, 0.0, 1.0, 0.0);
        testSingle(I, "f(x) = 1", [](Real x) -> Real { return 1.0; }, 0.0, 1.0, 1.0);
        testSingle(I, "f(x) = x", [](Real x) -> Real { return x; }, 0.0, 1.0, 0.5);
        testSingle(I, "f(x) = x^2",
                   [](Real x) -> Real { return x * x; }, 0.0, 1.0, 1.0/3.0);
        testSingle(I, "f(x) = sin(x)",
                   [](Real x) -> Real { return std::sin(x); }, 0.0, M_PI, 2.0);
        testSingle(I, "f(x) = cos(x)",
                   [](Real x) -> Real { return std::cos(x); }, 0.0, M_PI, 0.0);

        testSingle(I, "f(x) = Gaussian(x)",
                   NormalDistribution(), -10.0, 10.0, 1.0);
        testSingle(I, "f(x) = Abcd2(x)",
                   AbcdSquared(0.07, 0.07, 0.5, 0.1, 8.0, 10.0), 5.0, 6.0,
                   AbcdFunction(0.07, 0.07, 0.5, 0.1).covariance(5.0, 6.0, 8.0, 10.0));
    }

    template <class T>
    void testDegeneratedDomain(const T& I) {
        testSingle(I, "f(x) = 0 over [1, 1 + macheps]", [](Real x) -> Real { return 0.0; }, 1.0,
            1.0 + QL_EPSILON, 0.0);
    }

}


void IntegralTest::testSegment() {
    BOOST_TEST_MESSAGE("Testing segment integration...");

    using namespace integrals_test;

    testSeveral(SegmentIntegral(10000));
    testDegeneratedDomain(SegmentIntegral(10000));
}

void IntegralTest::testTrapezoid() {
    BOOST_TEST_MESSAGE("Testing trapezoid integration...");

    using namespace integrals_test;

    testSeveral(TrapezoidIntegral<Default>(integrals_test::tolerance, 10000));
    testDegeneratedDomain(TrapezoidIntegral<Default>(integrals_test::tolerance, 10000));
}

void IntegralTest::testMidPointTrapezoid() {
    BOOST_TEST_MESSAGE("Testing mid-point trapezoid integration...");

    using namespace integrals_test;

    testSeveral(TrapezoidIntegral<MidPoint>(integrals_test::tolerance, 10000));
    testDegeneratedDomain(TrapezoidIntegral<MidPoint>(integrals_test::tolerance, 10000));
}

void IntegralTest::testSimpson() {
    BOOST_TEST_MESSAGE("Testing Simpson integration...");

    using namespace integrals_test;

    testSeveral(SimpsonIntegral(integrals_test::tolerance, 10000));
    testDegeneratedDomain(SimpsonIntegral(integrals_test::tolerance, 10000));
}

void IntegralTest::testGaussKronrodAdaptive() {
    BOOST_TEST_MESSAGE("Testing adaptive Gauss-Kronrod integration...");

    using namespace integrals_test;

    Size maxEvaluations = 1000;
    testSeveral(GaussKronrodAdaptive(integrals_test::tolerance, maxEvaluations));
    testDegeneratedDomain(GaussKronrodAdaptive(integrals_test::tolerance, maxEvaluations));
}

void IntegralTest::testGaussLobatto() {
    BOOST_TEST_MESSAGE("Testing adaptive Gauss-Lobatto integration...");

    using namespace integrals_test;

    Size maxEvaluations = 1000;
    testSeveral(GaussLobattoIntegral(maxEvaluations, integrals_test::tolerance));
    // on degenerated domain [1,1+macheps] an exception is thrown
    // which is also ok, but not tested here
}

void IntegralTest::testTanhSinh() {
    BOOST_TEST_MESSAGE("Testing tanh-sinh integration...");

    using namespace integrals_test;
    testSeveral(TanhSinhIntegral());
}

void IntegralTest::testGaussLegendreIntegrator() {
    BOOST_TEST_MESSAGE("Testing Gauss-Legendre integrator...");

    using namespace integrals_test;

    const GaussLegendreIntegrator integrator(64);
    testSeveral(integrator);
    testDegeneratedDomain(integrator);
}

void IntegralTest::testGaussChebyshevIntegrator() {
    BOOST_TEST_MESSAGE("Testing Gauss-Chebyshev integrator...");

    using namespace integrals_test;

    const GaussChebyshevIntegrator integrator(64);
    testSingle(integrator, "f(x) = Gaussian(x)",
               NormalDistribution(), -10.0, 10.0, 1.0);
    testDegeneratedDomain(integrator);
}

void IntegralTest::testGaussChebyshev2ndIntegrator() {
    BOOST_TEST_MESSAGE("Testing Gauss-Chebyshev 2nd integrator...");

    using namespace integrals_test;

    const GaussChebyshev2ndIntegrator integrator(64);
    testSingle(integrator, "f(x) = Gaussian(x)",
               NormalDistribution(), -10.0, 10.0, 1.0);
    testDegeneratedDomain(integrator);
}



void IntegralTest::testGaussKronrodNonAdaptive() {
    BOOST_TEST_MESSAGE("Testing non-adaptive Gauss-Kronrod integration...");

    using namespace integrals_test;

    Real precision = integrals_test::tolerance;
    Size maxEvaluations = 100;
    Real relativeAccuracy = integrals_test::tolerance;
    GaussKronrodNonAdaptive gaussKronrodNonAdaptive(precision, maxEvaluations,
                                                    relativeAccuracy);
    testSeveral(gaussKronrodNonAdaptive);
    testDegeneratedDomain(gaussKronrodNonAdaptive);
}

void IntegralTest::testTwoDimensionalIntegration() {
    BOOST_TEST_MESSAGE("Testing two dimensional adaptive "
                       "Gauss-Lobatto integration...");

    using namespace integrals_test;

    const Size maxEvaluations = 1000;
    const Real calculated = TwoDimensionalIntegral(
        ext::shared_ptr<Integrator>(
            new TrapezoidIntegral<Default>(integrals_test::tolerance, maxEvaluations)),
        ext::shared_ptr<Integrator>(
            new TrapezoidIntegral<Default>(integrals_test::tolerance, maxEvaluations)))(
        std::multiplies<>(),
        std::make_pair(0.0, 0.0), std::make_pair(1.0, 2.0));

    const Real expected = 1.0;
    if (std::fabs(calculated-expected) > integrals_test::tolerance) {
        BOOST_FAIL(std::setprecision(10)
                   << "two dimensional integration: "
                   << "\n    calculated: " << calculated
                   << "\n    expected:   " << expected);
    }
}

namespace integrals_test {

    class sineF {
      public:
        Real operator()(Real x) const {
            return std::exp(-0.5*(x - M_PI_2/100));
        }
    };

    class cosineF {
      public:
        Real operator()(Real x) const {
            return std::exp(-0.5*x);
        }
    };

}

void IntegralTest::testFolinIntegration() {
    BOOST_TEST_MESSAGE("Testing Folin's integral formulae...");

    using namespace integrals_test;

    // Examples taken from
    // http://www.tat.physik.uni-tuebingen.de/~kokkotas/Teaching/Num_Methods_files/Comp_Phys5.pdf
    const Size nr[] = { 4, 8, 16, 128, 256, 1024, 2048 };
    const Real expected[] = { 4.55229440e-5,4.72338540e-5, 4.72338540e-5,
                              4.78308678e-5,4.78404787e-5, 4.78381120e-5,
                              4.78381084e-5};

    const Real t = 100;
    const Real o = M_PI_2/t;

    const Real tol = 1e-12;

    for (Size i=0; i < LENGTH(nr); ++i) {
        const Size n = nr[i];
        const Real calculatedCosine
            = FilonIntegral(FilonIntegral::Cosine, t, n)(cosineF(),0,2*M_PI);
        const Real calculatedSine
            = FilonIntegral(FilonIntegral::Sine, t, n)
                (sineF(), o,2*M_PI + o);

        if (std::fabs(calculatedCosine-expected[i]) > tol) {
            BOOST_FAIL(std::setprecision(10)
                << "Filon Cosine integration failed: "
                << "\n    calculated: " << calculatedCosine
                << "\n    expected:   " << expected[i]);
        }
        if (std::fabs(calculatedSine-expected[i]) > tol) {
            BOOST_FAIL(std::setprecision(10)
                << "Filon Sine integration failed: "
                << "\n    calculated: " << calculatedCosine
                << "\n    expected:   " << expected[i]);
        }
    }
}

namespace integrals_test {

    Real f1(Real x) {
        return 1.2*x*x+3.2*x+3.1;
    }

    Real f2(Real x) {
        return 4.3*(x-2.34)*(x-2.34)-6.2*(x-2.34) + f1(2.34);
    }

}

void IntegralTest::testDiscreteIntegrals() {
    BOOST_TEST_MESSAGE("Testing discrete integral formulae...");

    using namespace integrals_test;

    Array x(6), f(6);
    x[0] = 1.0; x[1] = 2.02; x[2] = 2.34; x[3] = 3.3; x[4] = 4.2; x[5] = 4.6;

    std::transform(x.begin(), x.begin()+3, f.begin(),   f1);
    std::transform(x.begin()+3, x.end(),   f.begin()+3, f2);

    const Real expectedSimpson =
        16.0401216 + 30.4137528 + 0.2*f2(4.2) + 0.2*f2(4.6);
    const Real expectedTrapezoid =
          0.5*(f1(1.0)  + f1(2.02))*1.02
        + 0.5*(f1(2.02) + f1(2.34))*0.32
        + 0.5*(f2(2.34) + f2(3.3) )*0.96
        + 0.5*(f2(3.3)  + f2(4.2) )*0.9
        + 0.5*(f2(4.2)  + f2(4.6) )*0.4;

    const Real calculatedSimpson =  DiscreteSimpsonIntegral()(x, f);
    const Real calculatedTrapezoid = DiscreteTrapezoidIntegral()(x, f);

    const Real tol = 1e-12;
    if (std::fabs(calculatedSimpson-expectedSimpson) > tol) {
        BOOST_FAIL(std::setprecision(16)
            << "discrete Simpson integration failed: "
            << "\n    calculated: " << calculatedSimpson
            << "\n    expected:   " << expectedSimpson);
    }

    if (std::fabs(calculatedTrapezoid-expectedTrapezoid) > tol) {
        BOOST_FAIL(std::setprecision(16)
            << "discrete Trapezoid integration failed: "
            << "\n    calculated: " << calculatedTrapezoid
            << "\n    expected:   " << expectedTrapezoid);
    }
}

void IntegralTest::testDiscreteIntegrator() {
    BOOST_TEST_MESSAGE("Testing discrete integrator formulae...");

    using namespace integrals_test;

    testSeveral(DiscreteSimpsonIntegrator(300));
    testSeveral(DiscreteTrapezoidIntegrator(3000));
}

namespace integrals_test{

std::vector<Real> x, y;

Real pw_fct(const Real t) { return QL_PIECEWISE_FUNCTION(x, y, t); }

void pw_check(const Integrator &in, const Real a, const Real b,
              const Real expected) {
    Real calculated = in(pw_fct, a, b);
    if (!close(calculated, expected))
        BOOST_FAIL(std::setprecision(16)
                   << "piecewise integration over [" << a << "," << b
                   << "] failed: "
                   << "\n   calculated: " << calculated
                   << "\n   expected:   " << expected
                   << "\n   difference: " << (calculated - expected));
}
} // empty namespace

void IntegralTest::testPiecewiseIntegral() {
    BOOST_TEST_MESSAGE("Testing piecewise integral...");

    using namespace integrals_test;

    x = { 1.0, 2.0, 3.0, 4.0, 5.0 };
    y = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
    ext::shared_ptr<Integrator> segment =
        ext::make_shared<SegmentIntegral>(1);
    ext::shared_ptr<Integrator> piecewise =
        ext::make_shared<PiecewiseIntegral>(segment, x);
    pw_check(*piecewise, -1.0, 0.0, 1.0);
    pw_check(*piecewise, 0.0, 1.0, 1.0);
    pw_check(*piecewise, 0.0, 1.5, 2.0);
    pw_check(*piecewise, 0.0, 2.0, 3.0);
    pw_check(*piecewise, 0.0, 2.5, 4.5);
    pw_check(*piecewise, 0.0, 3.0, 6.0);
    pw_check(*piecewise, 0.0, 4.0, 10.0);
    pw_check(*piecewise, 0.0, 5.0, 15.0);
    pw_check(*piecewise, 0.0, 6.0, 21.0);
    pw_check(*piecewise, 0.0, 7.0, 27.0);
    pw_check(*piecewise, 3.5, 4.5, 4.5);
    pw_check(*piecewise, 5.0, 10.0, 30.0);
    pw_check(*piecewise, 9.0, 10.0, 6.0);
}

namespace integrals_test {
    template <class T>
    void reportSiCiFail(
            const std::string& name, T z, T c, T e, Real diff, Real tol) {
        BOOST_FAIL(std::setprecision(16)
            << name << " calculation failed for " << z
            << "\n calculated: " << c
            << "\n expected:   " << e
            << "\n difference: " << diff
            << "\n tolerance:  " << tol);
    }
}

void IntegralTest::testExponentialIntegral() {
    BOOST_TEST_MESSAGE("Testing exponential integrals...");

    using namespace ExponentialIntegral;

    // reference values are calculated with Mathematica or Python/mpmath
    const Real data[][10] = {
            {0.0010000000000000000208, 0.0, 0.00099999994444444613193, 0.0, -6.330539864080593754, 0.0, -6.3295393640250381967, 0.0, 6.3315393641361493112, 0.0},
            {0.00070710678118654756077, 0.00070710678118654751747, 0.00070710682047025644818, 0.00070710674190283627304, -6.3305396140806145873, 0.785397913397448279, -6.329832507338711751, 0.7861055202179185354, 6.3312467208225174236, -0.78469130657697802259},
            {6.1232339957367660136e-20, 0.0010000000000000000208, 6.1232350162201617204e-20, 0.001000000055555557243, -6.330539364080593754, 1.570796326794896558, -6.3305398640805937539, 1.5717963267393410041, 6.330539864080593754, -1.5697963268504521119},
            {-0.00070710678118654747417, 0.00070710678118654760407, -0.00070710682047025636158, 0.00070710674190283635964, -6.3305396140806145873, 2.356194740192344837, -6.3312467208225174235, 2.3569013470128150935, 6.3298325073387117511, -2.3554871333718745805},
            {-0.0010000000000000000208, 1.2246467991473532027e-19, -0.00099999994444444613193, 1.2246465950321797194e-19, -6.330539864080593754, 3.141592653589793116, -6.3315393641361493112, 3.1415926535897931161, 6.3295393640250381967, -3.1415926535897931159},
            {-0.00070710678118654764736, -0.00070710678118654743088, -0.00070710682047025653477, -0.00070710674190283618645, -6.3305396140806145873, -2.3561947401923450819, -6.3312467208225174237, -2.3569013470128153382, 6.3298325073387117509, 2.3554871333718748256},
            {-1.8369701987210298041e-19, -0.0010000000000000000208, -1.8369705049015472569e-19, -0.001000000055555557243, -6.330539364080593754, -1.5707963267948968029, -6.3305398640805937541, -1.5717963267393412491, 6.3305398640805937538, 1.5697963268504523568},
            {0.00070710678118654738758, -0.00070710678118654769066, 0.00070710682047025627499, -0.00070710674190283644623, -6.3305396140806145873, -0.78539791339744852393, -6.3298325073387117512, -0.78610552021791878051, 6.3312467208225174234, 0.78469130657697826735},
            {0.10000000000000000555, 0.0, 0.099944461108276955702, 0.0, -1.7278683866572965838, 0.0, -1.6228128139692766136, 0.0, 1.8229239584193906159, 0.0},
            {0.07071067811865475853, 0.0707106781186547542, 0.070749950041603603669, 0.070671382625480302524, -1.7253704697591484327, 0.78289816362892975745, -1.6546990871336681256, 0.8586481132075703999, 1.7960418523846287393, -0.7171482131243632012},
            {6.123233995736766226e-18, 0.10000000000000000555, 6.1334444895968712791e-18, 0.10005557222505700112, -1.7228683861943336147, 1.5707963267948965577, -1.7278683866572965777, 1.670740787903173514, 1.7278683866572965899, -1.4708518656866196026},
            {-0.070710678118654749871, 0.07071067811865476286, -0.070749950041603595024, 0.070671382625480311198, -1.7253704697591484321, 2.3586944899608633586, -1.7960418523846287312, 2.4244444404654299234, 1.6546990871336681349, -2.2829445403822227075},
            {-0.10000000000000000555, 1.2246467991473532452e-17, -0.099944461108276955702, 1.2226067414369268591e-17, -1.7278683866572965838, 3.1415926535897931166, -1.8229239584193906159, 3.1415926535897931277, 1.6228128139692766136, -3.1415926535897931031},
            {-0.07071067811865476719, -0.070710678118654745541, -0.070749950041603612314, -0.07067138262548029385, -1.7253704697591484333, -2.3586944899608636035, -1.7960418523846287474, -2.4244444404654301511, 1.6546990871336681163, 2.2829445403822229697},
            {-1.8369701987210298678e-17, -0.10000000000000000555, -1.8400333469093536425e-17, -0.10005557222505700112, -1.7228683861943336147, -1.5707963267948968038, -1.7278683866572966021, -1.6707407879031737577, 1.7278683866572965654, 1.4708518656866198463},
            {0.070710678118654741211, -0.07071067811865477152, 0.070749950041603586379, -0.070671382625480319872, -1.7253704697591484315, -0.78289816362893000238, -1.6546990871336681441, -0.85864811320757066212, 1.7960418523846287232, 0.71714821312436342884},
            {0.25, 0.0, 0.2491335703197571641, 0.0, -0.82466306258094565309, 0.0, -0.54254326466191372953, 0.0, 1.0442826344437381945, 0.0},
            {0.17677669529663688651, 0.17677669529663687569, 0.17738935115398995617, 0.17616173766104897661, -0.80911938627521916259, 0.76977321991145556311, -0.63295764861417017313, 0.97841245803743094036, 0.98528112393626814822, -0.62363375572945104943},
            {1.5308084989341914715e-17, 0.25, 1.5468043259937507815e-17, 0.2508696848909122325, -0.79341294955282596203, 1.5707963267948965561, -0.82466306258094563794, 1.819929897114653724, 0.82466306258094566823, -1.3216627564751393958},
            {-0.17677669529663686486, 0.17677669529663689734, -0.17738935115398993475, 0.17616173766104899848, -0.80911938627521915876, 2.3718194336783375529, -0.98528112393626813018, 2.517958897860342088, 0.63295764861417019883, -2.1631801955523621542},
            {-0.25, 3.0616169978683829431e-17, -0.2491335703197571641, 3.0298246679137807606e-17, -0.82466306258094565309, 3.1415926535897931198, -1.0442826344437381945, 3.1415926535897931431, 0.54254326466191372953, -3.1415926535897930812},
            {-0.17677669529663690816, -0.17677669529663685404, -0.1773893511539899776, -0.17616173766104895473, -0.80911938627521916642, -2.3718194336783377978, -0.98528112393626816627, -2.5179588978603422901, 0.63295764861417014743, 2.163180195552362442},
            {-4.5924254968025744146e-17, -0.25, -4.6404129781529084775e-17, -0.2508696848909122325, -0.79341294955282596203, -1.5707963267948968087, -0.82466306258094569853, -1.8199298971146539613, 0.82466306258094560764, 1.3216627564751396331},
            {0.17677669529663684321, -0.17677669529663691899, 0.17738935115398991333, -0.17616173766104902036, -0.80911938627521915494, -0.769773219911455808, -0.63295764861417022453, -0.97841245803743122809, 0.98528112393626811213, 0.62363375572945125148},
            {0.5, 0.0, 0.49310741804306668916, 0.0, -0.17778407880661290134, 0.0, 0.45421990486317357992, 0.0, 0.55977359477616081175, 0.0},
            {0.35355339059327377302, 0.35355339059327375138, 0.3584268697133189464, 0.34860625536259795411, -0.11658254521497154353, 0.72290178026868503268, 0.23202371014762644078, 1.2063214162395304511, 0.46518880057756951253, -0.48946767681289259981},
            {3.0616169978683829431e-17, 0.5, 3.190788489546069124e-17, 0.50699674981966719583, -0.052776844956493615913, 1.5707963267948965502, -0.17778407880661287198, 2.0639037448379632547, 0.17778407880661293069, -1.0776889087518298763},
            {-0.35355339059327372973, 0.35355339059327379467, -0.35842686971331890493, 0.34860625536259799919, -0.11658254521497152822, 2.4186908733211080836, -0.46518880057756948275, 2.652124976776900558, -0.2320237101476263804, -1.9352712373502626237},
            {-0.5, 6.1232339957367658861e-17, -0.49310741804306668916, 5.8712695126944314728e-17, -0.17778407880661290134, 3.141592653589793131, -0.55977359477616081175, 3.1415926535897931642, -0.45421990486317357992, -3.1415926535897930366},
            {-0.35355339059327381632, -0.35355339059327370808, -0.35842686971331898787, -0.34860625536259790903, -0.11658254521497155883, -2.4186908733211083279, -0.4651888005775695423, -2.6521249767769007193, -0.23202371014762650116, 1.9352712373502629509},
            {-9.1848509936051488292e-17, -0.5, -9.5723654690421041556e-17, -0.50699674981966719583, -0.052776844956493615913, -1.5707963267948968264, -0.1777840788066129894, -2.0639037448379634696, 0.17778407880661281327, 1.0776889087518300913},
            {0.35355339059327368643, -0.35355339059327383797, 0.35842686971331886346, -0.34860625536259804427, -0.11658254521497151291, -0.72290178026868527697, 0.23202371014762632001, -1.2063214162395307784, 0.46518880057756945298, 0.48946767681289276116},
            {1.0, 0.0, 0.94608307036718301494, 0.0, 0.33740392290096813466, 0.0, 1.8951178163559367555, 0.0, 0.21938393439552027368, 0.0},
            {0.70710678118654754605, 0.70710678118654750275, 0.74519215535365930209, 0.66666481741950631398, 0.56680209825930888934, 0.53562961732242986806, 1.233466915678815284, 1.7803588648261259588, 0.099862719160197444251, -0.28997455411880742612},
            {6.1232339957367658861e-17, 1.0, 7.1960319005373041642e-17, 1.0572508753757285146, 0.83786694098020824089, 1.5707963267948965247, 0.33740392290096818619, 2.5168793971620796011, -0.33740392290096808314, -0.62471325642771357121},
            {-0.70710678118654745945, 0.70710678118654758935, -0.74519215535365923063, 0.66666481741950641427, 0.5668020982593089504, 2.605963036267363253, -0.099862719160197405024, 2.8516180994709857664, -1.2334669156788151226, -1.3612337887636670908},
            {-1.0, 1.2246467991473531772e-16, -0.94608307036718301494, 1.0305047480945999774e-16, 0.33740392290096813466, 3.1415926535897931723, -0.21938393439552027368, 3.1415926535897931934, -1.8951178163559367555, -3.1415926535897929056},
            {-0.70710678118654763265, -0.70710678118654741616, -0.74519215535365937355, -0.66666481741950621369, 0.56680209825930882828, -2.6059630362673634878, -0.099862719160197483478, -2.8516180994709858582, -1.2334669156788154453, 1.3612337887636674684},
            {-1.8369701987210297658e-16, -1.0, -2.158809570266204413e-16, -1.0572508753757285146, 0.83786694098020824089, -1.5707963267948969027, 0.33740392290096798009, -2.5168793971620797334, -0.33740392290096828924, 0.62471325642771370354},
            {0.70710678118654737286, -0.70710678118654767594, 0.74519215535365915917, -0.66666481741950651456, 0.56680209825930901147, -0.53562961732243010279, 1.2334669156788149613, -1.7803588648261263365, 0.099862719160197365796, 0.28997455411880751794},
            {1.5, 0.0, 1.3246835311721196804, 0.0, 0.47035631719539988668, 0.0, 3.301285449129797838, 0.0, 0.1000195824066326519, 0.0},
            {1.0606601717798213191, 1.0606601717798212541, 1.1839593855572855609, 0.9194789610047786165, 0.93002583013377829907, 0.22553329329273497242, 1.8495047911385570699, 2.5292224190594471214, -0.010546869128999664075, -0.16130364794487607552},
            {9.1848509936051488292e-17, 1.5, 1.3038076345658281264e-16, 1.7006525157682152449, 1.600632933361582593, 1.5707963267948964752, 0.47035631719539994775, 2.8954798579670162953, -0.4703563171953998256, -0.24611279562277693453},
            {-1.0606601717798211892, 1.060660171779821384, -1.1839593855572854849, 0.91947896100477878934, 0.93002583013377843491, 2.9160593602970581693, 0.010546869128999701077, 2.9802890056449171422, -1.8495047911385567612, -0.61237023453034594437},
            {-1.5, 1.8369701987210297658e-16, -1.3246835311721196804, 1.2215790424776667532e-16, 0.47035631719539988668, 3.1415926535897932298, -0.1000195824066326519, 3.1415926535897932111, -3.301285449129797838, -3.1415926535897926896},
            {-1.060660171779821449, -1.0606601717798211242, -1.1839593855572856369, -0.91947896100477844366, 0.93002583013377816324, -2.9160593602970583627, 0.010546869128999627072, -2.9802890056449171836, -1.8495047911385573786, 0.6123702345303462898},
            {-2.7554552980815446488e-16, -1.5, -3.9114229038065365107e-16, -1.7006525157682152449, 1.600632933361582593, -1.5707963267948970514, 0.47035631719539970344, -2.8954798579670163126, -0.47035631719540006991, 0.24611279562277695186},
            {1.0606601717798210593, -1.0606601717798215139, 1.1839593855572854089, -0.91947896100477896218, 0.93002583013377857075, -0.22553329329273516584, 1.8495047911385564526, -2.5292224190594474668, -0.010546869128999738079, 0.16130364794487611693},
            {2.0, 0.0, 1.6054129768026948486, 0.0, 0.4229808287748649957, 0.0, 4.9542343560018901634, 0.0, 0.048900510708061119567, 0.0},
            {1.4142135623730950921, 1.4142135623730950055, 1.6883194933582990243, 1.0649044719839209475, 1.1044891171909028954, -0.19981522706084708457, 2.1693935891748240916, 3.4589310472140426889, -0.039584645206981933184, -0.08229206049744467714},
            {1.2246467991473531772e-16, 2.0, 2.2208114946641807464e-16, 2.5015674333549756415, 2.4526669226469145219, 1.5707963267948963889, 0.42298082877486505138, 3.1762093035975914933, -0.42298082877486494002, 0.034616650007798203864},
            {-1.4142135623730949189, 1.4142135623730951787, -1.6883194933582989874, 1.0649044719839212109, 1.1044891171909031294, 3.3414078806506402814, 0.039584645206981962593, 3.0593005930923485567, -2.169393589174823594, 0.31733839362424952895},
            {-2.0, 2.4492935982947063545e-16, -1.6054129768026948486, 1.1135681831696612616e-16, 0.4229808287748649957, 3.1415926535897932894, -0.048900510708061119567, 3.1415926535897932219, -4.9542343560018901634, -3.1415926535897923336},
            {-1.4142135623730952653, -1.4142135623730948323, -1.6883194933582990613, -1.064904471983920684, 1.1044891171909026613, -3.3414078806506403646, 0.039584645206981903775, -3.059300593092348566, -2.1693935891748245892, -0.31733839362424937184},
            {-3.6739403974420595317e-16, -2.0, -6.6624344842268023737e-16, -2.5015674333549756415, 2.4526669226469145219, -1.5707963267948973103, 0.42298082877486482866, -3.1762093035975913914, -0.42298082877486516273, -0.03461665000779830579},
            {1.4142135623730947457, -1.4142135623730953519, 1.6883194933582989504, -1.0649044719839214744, 1.1044891171909033635, 0.19981522706084700137, 2.1693935891748230965, -3.458931047214042846, -0.039584645206981992002, 0.082292060497444686426},
            {5.0, 0.0, 1.5499312449446741373, 0.0, -0.19002974965664387862, 0.0, 40.185275355803177455, 0.0, 0.0011482955912753257973, 0.0},
            {3.5355339059327377302, 3.5355339059327375138, 3.6871508611543429709, -3.157181373909001357, -3.1547681046716106125, -2.1118502909279661892, -6.311949478580612776, 7.3697974788772077195, -0.0024132692373907452624, 0.0045042434314801607547},
            {3.0616169978683829431e-16, 5.0, 4.5436362172750887551e-15, 20.09321182569722639, 20.092063530105951065, 1.5707963267948920752, -0.19002974965664393734, 3.1207275717395707391, 0.1900297496566438199, -0.020865081850222464588},
            {-3.5355339059327372973, 3.5355339059327379467, -3.6871508611543449094, -3.1571813739090021642, -3.1547681046716114182, 5.2534429445177613695, 0.0024132692373907438925, 3.1460968970212734025, 6.3119494785806111631, 4.2282048252874106008},
            {-5.0, 6.1232339957367658861e-16, -1.5499312449446741373, -1.174343542809398959e-16, -0.19002974965664387862, 3.1415926535897932037, -0.0011482955912753257973, 3.1415926535897932376, -40.185275355803177455, -3.1415926535897750631},
            {-3.5355339059327381632, -3.5355339059327370808, -3.6871508611543410324, 3.1571813739090005499, -3.1547681046716098067, -5.2534429445177574859, 0.0024132692373907466323, -3.1460968970212733959, 6.3119494785806143889, -4.2282048252874183613},
            {-9.1848509936051488292e-16, -5.0, -1.3630908648516543815e-14, -20.09321182569722639, 20.092063530105951065, -1.5707963267949102514, -0.19002974965664370247, -3.1207275717395708086, 0.19002974965664405477, 0.020865081850222534065},
            {3.5355339059327368643, -3.5355339059327383797, 3.6871508611543468479, 3.1571813739090029713, -3.154768104671612224, 2.1118502909279700728, -6.3119494785806095501, -7.3697974788771999589, -0.0024132692373907425226, -0.004504243431480167346},
            {10.0, 0.0, 1.6583475942188740493, 0.0, -0.045456433004455372635, 0.0, 2492.2289762418777591, 0.0, 4.1569689296853242774e-6, 0.0},
            {7.0710678118654754605, 7.0710678118654750275, -3.7745175303432929615, 62.642575559232345941, 62.642571122904060317, 5.3452347019798827768, 125.28514668213645736, -7.5489559055283299711, 4.4363282856236060522e-6, -0.000079155158306803868409},
            {6.1232339957367658861e-16, 10.0, 6.7436601941373126554e-13, 1246.1144901994233444, 1246.1144860424544147, 1.5707963267942222532, -0.045456433004455405946, 3.2291439210137707199, 0.045456433004455339323, 0.087551267423977378721},
            {-7.0710678118654745945, 7.0710678118654758935, 3.7745175303433438135, 62.642575559232397046, 62.642571122904111422, -2.2036420483901403904, -4.4363282856235323221e-6, 3.1415134984314864345, -125.28514668213635515, -10.690548559118021506},
            {-10.0, 1.2246467991473531772e-15, -1.6583475942188740493, -6.6623371218859982286e-17, -0.045456433004455372635, 3.1415926535897933412, -4.1569689296853242774e-6, 3.1415926535897932385, -2492.2289762418777591, -3.1415926535870957744},
            {-7.0710678118654763265, -7.0710678118654741616, 3.7745175303432421091, -62.642575559232294835, 62.642571122904009212, 2.2036420483900386859, -4.4363282856236797822e-6, -3.1415134984314864347, -125.28514668213655957, 10.690548559118224914},
            {-1.8369701987210297658e-15, -10.0, -2.0230980582411937966e-12, -1246.1144901994233444, 1246.1144860424544147, -1.5707963267969197173, -0.045456433004455272699, -3.2291439210137705144, 0.04545643300445547257, -0.087551267423977584235},
            {7.0710678118654737286, -7.0710678118654767594, -3.7745175303433946658, -62.642575559232448152, 62.642571122904162528, -5.3452347019799844813, 125.28514668213625294, 7.5489559055281265623, 4.4363282856234585915e-6, 0.000079155158306804015139},
            {15.0, 0.0, 1.6181944437083687391, 0.0, 0.046278677674360439604, 0.0, 234955.85249076830358, 0.0, 1.9186278921478669771e-8, 0.0},
            {10.606601717798213191, 10.606601717798212541, -471.05515346873571713, -1328.2594907541181068, -1328.2594913017652105, 472.62595127241629118, -2656.5189820558856064, -942.11030841435617353, 5.4764710367805069957e-7, 1.4768856774331466472e-6},
            {9.1848509936051488292e-16, 15.0, 1.0008479153886857898e-10, 117477.92624539374493, 117477.92624537455865, 1.5707963266948118277, 0.046278677674360479423, 3.1889907705032654049, -0.046278677674360399786, 0.047398116913472073375},
            {-10.606601717798211892, 10.60660171779821384, 471.05515346873477896, -1328.2594907541203959, -1328.2594913017674995, -469.48435861882555978, -5.4764710367805350438e-7, 3.1415941304754706716, 2656.5189820558810283, -945.2519010679478431},
            {-15.0, 1.8369701987210297658e-15, -1.6181944437083687391, 7.9637292204322254463e-17, 0.046278677674360439604, 3.1415926535897933315, -1.9186278921478669771e-8, 3.1415926535897932385, -234955.85249076830358, -3.1415926531894540723},
            {-10.60660171779821449, -10.606601717798211242, 471.05515346873665531, 1328.2594907541158178, -1328.2594913017629215, 469.48435861882743613, -5.4764710367804789476e-7, -3.1415941304754706716, 2656.5189820558901845, 945.25190106794409045},
            {-2.7554552980815446488e-15, -15.0, -3.0025437461660077385e-10, -117477.92624539374493, 117477.92624537455865, -1.5707963270951509938, 0.046278677674360320148, -3.1889907705032652188, -0.04627867767436055906, -0.047398116913472259445},
            {10.606601717798210593, -10.606601717798215139, -471.05515346873384079, 1328.2594907541226849, -1328.2594913017697886, -472.62595127241441484, -2656.5189820558764503, 942.11030841435992621, 5.4764710367805630921e-7, -1.4768856774331489463e-6},
            {20.0, 0.0, 1.5482417010434398402, 0.0, 0.04441982084535331654, 0.0, 25615652.66405658882, 0.0, 9.8355252906498816904e-11, 0.0},
            {14.142135623730950921, 14.142135623730950055, 24486.68965358318626, 26235.092183360235647, 26235.092183384186671, -24485.118857281672279, 52470.184366744507202, 48973.379307191653857, -2.3951024244493564362e-8, -2.5280915398646225218e-8},
            {1.2246467991473531772e-15, 20.0, 1.4853900090407493933e-8, 12807826.332028294459, 12807826.332028294361, 1.5707963119409965288, 0.044419820845353372442, 3.1190380278383364344, -0.044419820845353260638, -0.02255462575145675408},
            {-14.142135623730949189, 14.142135623730951787, -24486.689653583186682, 26235.092183360320531, 26235.092183384271555, 24488.260449935262494, 2.3951024244493652701e-8, 3.1415926283088778398, -52470.184366744337435, 48970.23771453806322},
            {-20.0, 2.4492935982947063545e-15, -1.5482417010434398402, 1.1180354790034662799e-16, 0.04441982084535331654, 3.1415926535897931885, -9.8355252906498816904e-11, 3.1415926535897932385, -25615652.66405658882, -3.1415925941741928768},
            {-14.142135623730952653, -14.142135623730948323, -24486.689653583185838, -26235.092183360150763, 26235.092183384101787, -24488.26044993526165, 2.3951024244493476023e-8, -3.1415926283088778398, -52470.184366744676971, -48970.237714538064908},
            {-3.6739403974420595317e-15, -20.0, -4.4561700271222483452e-8, -12807826.332028294459, 12807826.332028294361, -1.5707963713565968905, 0.044419820845353148834, -3.1190380278383365344, -0.044419820845353484245, 0.022554625751456854031},
            {14.142135623730947457, -14.142135623730953519, 24486.689653583187103, -26235.092183360405415, 26235.09218338435644, 24485.118857281673122, 52470.184366744167666, -48973.37930719165217, -2.3951024244493741041e-8, 2.528091539864622434e-8},
            {24.989999999999998437, 0.0, 1.5315374843262580141, 0.0, -0.0072448862399538447391, 0.0, 2977286754.0962403117, 0.0, 5.4047414055481939457e-13, 0.0},
            {17.67059846185182207, 17.670598461851820988, -885673.16331345207799, -400248.00540652140574, -400248.00540652215797, 885674.7341097792085, -800496.01081304623667, -1771346.3266269055961, 7.5223231588352466724e-10, 3.3560562518531656391e-10},
            {1.5301961755346176992e-15, 24.989999999999998437, 2.1825789542468514196e-6, 1488643377.0481201559, 1488643377.0481201559, 1.5707941442159423724, -0.0072448862399538534498, 3.1023338111211545727, 0.0072448862399538360284, -0.039258842468638544566},
            {-17.670598461851819906, 17.670598461851823152, 885673.16331345318248, -400248.0054065240787, -400248.00540652483093, -885671.59251712672318, -7.5223231588352706351e-10, 3.1415926539253988636, 800496.01081304089075, -1771349.4682195569769},
            {-24.989999999999998437, 3.0603923510692353985e-15, -1.5315374843262580141, -1.7421457417638512695e-17, -0.0072448862399538447391, 3.1415926535897931172, -5.4047414055481939457e-13, 3.1415926535897932385, -2977286754.0962403117, -3.1415839232739762511},
            {-17.670598461851824234, -17.670598461851818824, 885673.16331345097349, 400248.00540651873277, -400248.005406519485, 885671.5925171245142, -7.5223231588352227096e-10, -3.1415926539253988636, 800496.01081305158263, 1771349.4682195613948},
            {-4.5905885266038530977e-15, -24.989999999999998437, -6.5477368627405542588e-6, -1488643377.0481201559, 1488643377.0481201559, -1.5708028745317593598, -0.0072448862399538186069, -3.1023338111211548151, 0.0072448862399538708713, 0.039258842468638787005},
            {17.670598461851817742, -17.670598461851825316, -885673.16331345428695, 400248.00540652675168, -400248.00540652750391, -885674.73410978141745, -800496.01081303554482, 1771346.3266269011781, 7.5223231588352945979e-10, -3.3560562518531458362e-10}
    };

    const Real tol = 1e-10;

    for (const auto& i : data) {
        const Real x = i[0];
        const Real y = (std::abs(i[1]) < 1e-12) ? 0.0 : i[1];
        const std::complex<Real> z(x, y);


        const std::complex<Real> si = Si(z);
        std::complex<Real> ref(i[2], i[3]);
        Real diff = std::abs(si-ref)/std::abs(ref);
        if (diff > tol
            || (std::abs(ref.real()) < tol && std::abs(si.real()) > tol)
            || (std::abs(ref.imag()) < tol && std::abs(si.imag()) > tol)) {
            integrals_test::reportSiCiFail("Si", z, si, ref, diff, tol);
        }

        const std::complex<Real> ci = Ci(z);
        ref = std::complex<Real>(i[4], i[5]);
        diff = std::min(std::abs(ci-ref), std::abs(ci-ref)/std::abs(ref));
        if (diff > tol
            || (std::abs(ref.real()) < tol && std::abs(ci.real()) > tol)
            || (std::abs(ref.imag()) < tol && std::abs(ci.imag()) > tol)) {
            integrals_test::reportSiCiFail("Ci", z, ci, ref, diff, tol);
        }

        const std::complex<Real> ei = Ei(z);
        ref = std::complex<Real>(i[6], i[7]);
        diff = std::abs(ei-ref)/std::abs(ref);
        if (diff > tol
            || (std::abs(ref.real()) < tol && std::abs(ei.real()) > tol)
            || (std::abs(ref.imag()) < tol && std::abs(ei.imag()) > tol)) {
            integrals_test::reportSiCiFail("Ei", z, ei, ref, diff, tol);
        }

        const std::complex<Real> e1 = E1(z);
        ref = std::complex<Real>(i[8], i[9]);
        diff = std::min(std::abs(e1-ref), std::abs(e1-ref)/std::abs(ref));
        if (std::abs(z) < 10.0)
            if (diff > tol
            || (std::abs(ref.real()) < tol && std::abs(e1.real()) > tol)
            || (std::abs(ref.imag()) < tol && std::abs(e1.imag()) > tol)) {
            integrals_test::reportSiCiFail("E1", z, e1, ref, diff, tol);
        }
    }
}



void IntegralTest::testRealSiCiIntegrals() {
    BOOST_TEST_MESSAGE("Testing real Ci and Si...");

    using namespace ExponentialIntegral;

    // reference values are calculated with Mathematica or Python/mpmath
    const Real data[][3] = {
            {1e-12, 1e-12, -27.0538054510270153677},
            {0.1, 0.09994446110827695570, -1.7278683866572965838},
            {1.0, 0.9460830703671830149, 0.3374039229009681347},
            {1.9999, 1.6053675097543679041, 0.4230016343635392},
            {3.9999, 1.758222058430840841, -0.140965355646150101},
            {4.0001, 1.758184218306157867, -0.140998037827177150},
            {5.0, 1.5499312449446741373, -0.19002974965664387862},
            {7.0, 1.4545966142480935906, 0.076695278482184518383,},
            {10.0, 1.6583475942188740493, -0.045456433004455372635},
            {15.0, 1.6181944437083687391, 0.046278677674360439604},
            {20.0, 1.5482417010434398402, 0.04441982084535331654},
            {24.9, 1.532210740207620024, -0.010788215638781789846},
            {25.1, 1.5311526281483412938, -0.0028719014454227088097},
            {30.0, 1.566756540030351111, -0.033032417282071143779},
            {40.0, 1.5869851193547845068, 0.019020007896208766962},
            {400.0, 1.5721148692738117518, -0.00212398883084634893},
            {4000.0, 1.5709788562309441985, -0.00017083030544201591130}
    };


    const Real tol = 1e-12;

    for (const auto& i : data) {
        Real x = i[0];
        Real si = Si(x);

        Real diff = std::fabs(si - i[1]);
        if (diff > tol) {
            integrals_test::reportSiCiFail("SineIntegral", x, si, i[1], diff, tol);
        }

        const Real ci = Ci(x);
        diff = std::fabs(ci - i[2]);
        if (diff > tol) {
            integrals_test::reportSiCiFail("CosineIntegral", x, ci, i[2], diff, tol);
        }

        x = -i[0];
        si = Si(x);
        diff = std::fabs(si + i[1]);
        if (diff > tol) {
            integrals_test::reportSiCiFail("SineIntegral", x, si, Real(-i[1]), diff, tol);
        }
    }
}

test_suite* IntegralTest::suite() {
    auto* suite = BOOST_TEST_SUITE("Integration tests");
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testSegment));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTrapezoid));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testMidPointTrapezoid));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testSimpson));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussKronrodAdaptive));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussKronrodNonAdaptive));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussLobatto));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussLegendreIntegrator));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussChebyshevIntegrator));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testGaussChebyshev2ndIntegrator));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTwoDimensionalIntegration));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testFolinIntegration));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testDiscreteIntegrals));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testDiscreteIntegrator));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testPiecewiseIntegral));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testExponentialIntegral));
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testRealSiCiIntegrals));

#ifdef QL_BOOST_HAS_TANH_SINH
    suite->add(QUANTLIB_TEST_CASE(&IntegralTest::testTanhSinh));
#endif

    return suite;
}

